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We present a new fully first order strongly hyperbolic representation of the BSSN formulation of 
Einstein's equations with optional constraint damping terms. We describe the characteristic fields 
of the system, discuss its hyperbolicity properties, and present two numerical implementations and 
simulations: one using finite differences, adaptive mesh refinement and in particular binary black 
holes, and another one using the discontinuous Galerkin method in spherical symmetry. The results 
of this paper constitute a first step in an effort to combine the robustness of BSSN evolutions with 
very high accuracy numerical techniques, such as spectral collocation multi-domain or discontinuous 
Galerkin methods. 
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I. INTRODUCTION 

Complete, long term numerical simulations of the in- 
spiral, merger and ringdown of two black holes became 
possible a few years ago [1-3] and are now carried out 
by numerous groups; see [4-6] for recent reviews on the 
topic. One of the motivations for studying the dynam- 
ics of these inspiraling compact binaries is due to the 
fact that they are among the most promising sources of 
gravitational waves for the upcoming advanced network 
of earth based laser interferometric detectors [7]. More- 
over, with (and only with) modeling of enough accuracy, 
these detectors should be able to extract from the waves 
physical data about these sources such as the component 
masses and spins. 

Until a few years ago, such simulations were plagued 
by short-term instabilities. With full, long and stable 
simulations now being carried out systematically, devel- 
opment efforts have focused on efficiency and accuracy, 
better boundary conditions (see [8] for a review), and 
wave extraction methods (see, for example, [9-12] and 
references therein), all of which are especially impor- 
tant for many-orbit evolutions, such as those needed to 
make comparisons with post-Newtonian models, and cal- 
ibration or fitting of semi-analytical or phenomcnologi- 
cal models [13-33]. We consider here only solutions of 
the vacuum Einstein equations, and intentionally ignore 
the much larger and astrophysically probably even more 



interesting case where matter, radiation, or electromag- 
netic fields are present. 

Most, if not all, numerical simulations of binary black 
holes (BBH) currently use one of two formulations of the 
Einstein equations. One of them is the generalized har- 
monic system, which has been successfully implemented 
in BBH simulations using finite difference adaptive mesh 
refinement (AMR) [4], pseudospectral collocation codes 
(see, for example, [34-38] and references therein), and 
multi-domain finite differences [39], in a first order in 
time and second order in space formulation [40] in the 
AMR case, and a fully first order reduction [41] other- 
wise. In either case the key ingredient is a constraint 
damping mechanism [42], originally proposed in [43] (in 
that reference referred to as A-systems because they were 
first introduced as Lagrange multipliers). The other one 
is the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) 
system [44, 45], which has been implemented by many 
groups using finite difference codes in a first order in 
time, second order in space form (see [5] for a review)- 
we refer to this as simply BSSN or second-order BSSN 
(as opposed to our fully first order reduction, to which 
we will refer as FOBSSN). Some variant of the "stan- 
dard gauge" conditions for the BSSN formulation, con- 
sisting of 1+log slicing or generalizations thereof and the 
so called Gamma-driver shift [46] condition are, more of- 
ten than not, used. The hyperbolicity of BSSN with a 
generalization of these gauge conditions was studied in 
Refs. [47, 48]. 
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While the Einstein equations are fundamentally a sec- 
ond order system, many advanced numerical techniques 
for hyperbolic systems - such as multi-domain high order 
finite difference, spectral collocation, and discontinuous 
Galerkin methods - are well developed for first order hy- 
perbolic systems. At the same time, the standard second 
order in space BSSN system with the standard gauge 
conditions has shown remarkably robust properties in a 
variety of compact binary configurations. One wonders, 
then, if it is possible to combine some of the numeri- 
cal techniques that are often used for very high accuracy 
simulations of hyperbolic differential equations with the 
BSSN system. 

One recent approach has been to adapt advanced tech- 
niques for fully first order systems to second order in 
space ones [49-53]. Perhaps paradoxically, it appears to 
be more difficult to guarantee stability for naturally sec- 
ond order systems than for first order reductions of them, 
though progress is being made on this front (see, for ex- 
ample, [52, 54, 55]). Another approach is to rewrite the 
Einstein equations as a fully first order hyperbolic sys- 
tem. In this paper we explore the latter and we refer to 
our first order reduction of BSSN as FOBSSN. 

This paper is organized as follows. Section II reviews 
the BSSN system in covariant form. The first order re- 
duction is carried out in Section III, where we also show 
that FOBSSN is strongly hyperbolic under suitable con- 
ditions on the gauge parameters, and discuss the prop- 
agation of the constraints. Section IV summarizes some 
of our results from numerical simulations of binary black 
holes using FOBSSN and adaptive-mesh refinement and 
finite differences, and a multi-domain nodal discontinu- 
ous Galerkin scheme in spherical symmetry. When ap- 
propriate, we compare results from our FOBSSN simu- 
lations to simulations using BSSN in its standard form. 
A preliminary look at turduckening [56-58] for a polyno- 
mial/spectral Galerkin method is presented in the con- 
text of the spherically reduced FOBSSN system. Appen- 
dices collect further details on the covariant BSSN system 
as well as expressions for the fundamental fields in terms 
of characteristic variables. 



II. REVIEW OF THE BSSN SYSTEM 

We briefly review the second order form of BSSN [44, 
45] with moving puncture gauge conditions. Here we fol- 
low the approach of Ref . [59] , which is spatially-covariant 
(but not fully space-time covariant). The spatial metric 
and extrinsic curvature are denoted 7^ and Kij, respec- 
tively. These are replaced by the BSSN variables 
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where 7 is the determinant of 7^^ and 7 is a fiducial scalar 
density of weight 2 which remains to be specified. Also, 
K — ^"^^Kij is the trace of the extrinsic curvature. The 
variable 7^^ is the conformal metric, (/> is the conformal 
factor, and is the conformally rescaled, trace-free 
part of the extrinsic curvature. (In addition to this "(/)" 
variant of BSSN, there exist also the W and % variants 
where the variable is replaced by = or 
X = (7/7)"^^"^! respectively.) These variables are re- 
stricted by the algebraic conditions 



7 = 7 







(2) 



where 7 is the determinant of 7ij.^ The BSSN system 
also includes the "conformal connection vector" 



A^ = 7J'=Ar 



(3) 



as independent variables, and we have defined AF^jfe = 

f%fc — r'jfc. Here, f 'j-fc are the Christoffel symbols built 

from the conformal metric and F^jfc is a fiducial connec- 
tion. In this covariant language the BSSN variables are 
tensors with no density weights. In particular, </) is a 
scalar and A' is a contravariant vector. 

It is. often convenient to consider the fiducial connec- 
tion F jfe to be constructed from a "fiducial metric" 7^^ 
whose determinant is 7. We stress that the fiducial fields 
are not dynamical variables. They are freely chosen func- 
tions, required by covariance. Throughout the main body 
of the paper we assume that the fiducial connection is 
built from a flat, time independent metric 7^^ whose 
determinant is 7. If the coordinates are interpreted as 
Cartesian, then 7^^ = diag(l,l,l). In this case 7=1 



and F 



0. The vector A' then reduces to the con- 



formal connection functions F* = T"^ jkl^^ and Eqs. (4) 
below reduce to the usual second order BSSN system. 
The evolution equations for the BSSN variables are 
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^ Note that we use both 7 and 7 in our notation. 
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where a is the lapse function and is the shift vec- 
tor. Also, the time derivative operator is defined by 
= dt — Cfj, where is the Lie derivative with re- 
spect to the shift. Next, Di, Di and Di are the covari- 
ant derivatives built from the physical metric, conformal 
metric, and fiducial metric, respectively, and DiDja — 
'Di'Dja - AT''ijDka. Finally, in Eq. (4d), TF denotes 
the trace-free part of the expression in brackets. 

The Ricci tensor can be written as a sum of two terms, 

R,^ = R^J + R% ■ (5) 
The Ricci tensor for the conformal metric is 

R^J = -^T^Sfc^n^j +7fc(«^j)^'''+7'"Af^„Af(y)fc 

+7'=^[2Af'"fe(,Af,)„, + Af™,fcAf„,,] , (6) 

and the term i?^ is defined by 

Rf^ - -2D,Dj(j)~2%D''Dk(l) 

+AD,4>D,cb - A%&(l,Di^ . (7) 

All tensors with a tilde have their indices raised and low- 
ered with the conformal metric7ij . Details of the deriva- 
tion of the equations of motion (4) are contained in Ap- 
pendix A. 

In addition to the algebraic constraints, solutions to 
the second order BSSN system must satisfy a set of differ- 
ential constraints stemming from the 3-1-1 decomposition 
(Hamiltonian and momentum constraints) and definition 
of the conformal connection vector. Expressed in terms 
of the evolved variables, these are given by 

n = e-^'>' (^R - 8D'D,(t) - 8D'(l}D,(j)^ 

M' = DjA'^ + 6A'^dj(l) - \]D'K = , (8b) 

^ l'^Af]t, = , (8c) 

where R = 'j'^^Rij. Using the Bianchi identities and the 
BSSN equations (4), one can derive a closed homoge- 
neous evolution system for the constraint fields H, 
and g^. This constraint evolution system can be written 
in first order symmetric hyperbolic form [47, 58, 60, 61]. 
Therefore, if the initial data satisfies the constraints, then 
the constraints will be preserved for all times as long as 
suitable boundary conditions are provided. Constraint- 
preserving boundary conditions for the BSSN system 
have been discussed in Refs. [60, 61]. 

Black hole evolutions with the second-order BSSN sys- 
tem are typically carried out using the moving puncture 
gauge conditions consisting of H-log slicing and Gamma- 
driver shift. In this paper we consider the general Bona- 
Masso slicing condition [62], written in the form [47] 

dj_a = -aV(a, <I^)K + So, , (9) 
where f{a,(f)) is an arbitrary positive function of the lapse 
a and the conformal factor (j). The source term Sa is a 



function of the spatial coordinates. For 1 -I- log slicing, 
/(a,(/)) = 2/a and = 0. 

The shift condition considered in this paper is a gen- 
eralization of the Gamma-driver shift, written as [47] 

do(3' = a^G{a, <j>)B' + S'p , (10a) 
doB' = e-4^i?(a, 4,)doK' - 7]B' + S'b . (10b) 

Here, is an auxiliary field and the time derivative op- 
erator is defined by do = dt — fS^Dj. The functions G and 
H depend on the lapse a and conformal factor (p. The 
source terms S*^ and Sg are functions of the spatial co- 
ordinates. The standard choices for the Gamma-driver 
shift condition are G(a,0) = 3/(4q!2), i/(a,0) = e'*'^, 
Sj^ = S'b = 0, and 7] = 3/(4M), where M is the ADM 
mass of the system or another relevant mass scale. (If dif- 
ferent regions of the domain have different mass scales, 
e.g. for binary black hole systems with a large mass ratio, 
then r/ may vary in space [63, 64].) 

III. FIRST-ORDER BSSN 

A. Evolution System 

The BSSN system as described above contains second 
order derivatives in space acting on the variables a, /3% 
(j) and "fij. To write the system in fully first order form 
we introduce the new variables 

= Dia , (11a) 

A^" = , (lib) 

0, = D,^ , (11c) 

%i] = Dk%j ■ (lid) 

These definitions yield the associated constraints 

Ai=ai- Dia = , (12a) 

= - D,p^ = , (12b) 

= 0, - = , (12c) 

T^kij = Iki] - Dk%j = . (12d) 

Observe that the derivative Dk applied to the algebraic 
constraint 7 = 7 yields the condition 

f'lk^J = . (13) 

This is a new algebraic constraint that the first order 
BSSN variables must satisfy, along with the algebraic 
constraints (2) inherited from second order BSSN. 

The evolution equations for the new variables are ob- 
tained by computing their time derivatives using the 
second-order BSSN equations (4a), (4b) and gauge con- 
ditions (9), (10a). In carrying out these calculations we 
continue to assume that the fiducial metric is flat and 
time-independent. The complete system of first order 
equations, the FOBSSN system, can be conveniently split 
into gauge and non-gauge sectors. The gauge sector is 
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doa = 

doB^ = 
W = 



-a^fm - 2faa,K - a" [f^a, + f^^,] K + + - ^"A, , 

a^GB^ + Si , 

a^GD.B^ + 2GaaiB^ + [Gm + G^<l>^] B^ + f3,''/3k^ + D.Sj; - k^^B^^ 

1 



(14a) 
(14b) 
(14c) 
(14d) 
(14e) 

(14f) 



r 



Here, subscripts a and (p on the functions / and G denote 
partial derivatives. Note that terms proportional to the 
constraints Ai and Bi^ have been added to the right- 
hand sides of the evolution equations for and /3i^ . The 
corresponding proportionality constants are and . 



These terms can be used as a damping mechanism for 
any numerical violation of the constraints Ai — Q and 
B,^ = 0. 

The remaining evolution equations, which comprise the 
non-gauge sector, are 



da A, J 



6 6 6 

I 2 I. 
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+aKA,j - 2aA,kA) + 2 Am (3. 



ijPk 



(15a) 
(15b) 
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Here, we have defined 



-2aDkAj + 2(Dfe/3(/)7,), - -j,jDk(3/ 

2 

-2akAij + /3k lUj + 2^ki(iPj) - ^ Tfcij/?^ - K?T^kij , 
-AT'^kl'^Pi' + ^Ar.kf'P/ - 2A'^aj + 2a (APfe^i" + 6i^^0,) 



(15e) 



(15f) 



AT\i = -7'^ [^ki] + Ilk] - l]ki) , 



R.] = --r'^fc7^.,+7M«^j)A'^ + f"ArLAr(,,), + 7'=^[2Ar™,(,Ar,),„, + Ar™ , 



which follow from the definition of the Christoffel sym- 
bols and the identity (6) for the Ricci tensor. In the 
evolution equations for (/>i and "fkiji the constraints Ci 
and T>kij are subtracted with constants k'^ and k"^ . These 
terms are included as damping mechanisms for these con- 
straints, in analogy with k" and . The term propor- 
tional to the constant a in the evolution equation for 



A\ Eq. (15f), is equal to the constraint 2^^^D^Bk]'' = 
0. This term is needed to make the evolution system 
strongly hyperbolic, as discussed below. 
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B. Constraint Propagation 

The FOBSSN system is subject to the algebraic con- 
straints 7-7 = 0, A] = f^Aij = 0, and 7^/ = ^'^^^j = 
0. As discussed in the next section, our numerical codes 
enforce some but not all of these constraints. If the al- 
gebraic constraints are not enforced, but free to evolve, 
the first order BSSN evolution equations (15c), (15d) and 
(15c) imply 



do ln(7/7) 
doAl 



-2aAl , (16a) 
aKAl , (16b) 

~2aDkAl - 2akAl + h'lu' 



(16c) 



It follows from Eq. (16b) that A\ is zero along an inte- 
gral curve c oi do = dt — 13^ Dj if it is zero at some point 
on this curve. Therefore, if all integral curves c intersect 
the initial surface, it is sufficient to require A* = on 
this surface in order to guarantee that the algebraic con- 
straint A\ = Q holds at every time and everywhere on the 
computational domain. On the other hand, if c intersects 
the boundary surface, which happens if the shift vector 
is pointing outwards at the boundary, then A\ = needs 
to be enforced as a boundary condition in order to guar- 
antee the satisfaction of the algebraic constraint A\ = 0. 
If this constraint holds, it follows by a similar argument 
from Eq. (16a) that the determinant constraint 7 = 7 
holds if it is satisfied initially and suitable boundary con- 
ditions are specified in case (3^ is outward pointing at the 
boundary. Eqs. (16a, 16b) also show that is is consistent 
to enforce the algebraic constraints 7 = 7 and A\ = Q 
throughout evolution, as is the case for the second order 
BSSN system. 

On the other hand, it does not follow immediately from 



Eq. (16c) and suitable initial and boundary conditions 
that the trace constraint 7^^' = holds, unless the term 
A'Wkij is zero.Notice that ^'^Vkij = ^ki" - ^^fc ln(7/7) 
so this term can be expressed in terms of the algebraic 
constraints. This means that the propagation of the al- 
gebraic trace constraint 7^^' = is coupled to those of 
the constraints Ai = 0, Bi^ = 0, Q = and Vkij = 0, 
and one cannot consistently enforce 7/ci* = along with 
the other algebraic constraints unless A'^^Vkij = 0. 

Alternatively, it is possible to decouple the algebraic 
constraints from the remaining ones by adding the term 

- —j^^A'^'Vklm + y7.,f (17) 

to the right-hand side of Eq. (15c), which has the same 
effect as the replacements 



1 



DkA^, ^ D,A, - 



and 



+ 37.. 



in that equation. With this, the last two terms on the 
right-hand side of Eq. (16c) drop, and one obtains a 
closed, homogeneous evolution system for the algebraic 
constraints. Therefore, it is consistent to set the alge- 
braic constraints to zero even if 'Dkij 7^ 0- 

We now consider the constraints Ai = 0, Bi^ = 0, Ci — 
and Vkij = that were introduced in the reduction 
of BSSN to first-order. The evolution equations imply 
that the constraint fields Ai , Bi^ , and Vkij satisfy the 
following linear, homogeneous system of equations: 



doA^ = -2afKA^ - a^K + + {D,P^)Aj + a^B^ - k"A^ , (18a) 

doBi = 2aGB'A^ + a^B' [G^A^ + G^Q + (D,^'')Bk' + fik'B,'' - , (18b) 

doC, = -^A + {D,(3^)Cj + c^jB^ - K-^C, , (18c) 

^oDk^, = ~2A,jAk + (Pk[3^)Vuj + lUjBk^ + 2I?fe£(,/3,)^ - ^ Phj/?/ - ^^^Vk^J . (18d) 



(The term in Eq. (17) has to be added to the right-hand 
side of the last equation in case Eq. (15e) is modified in 
the way described above.) 

Here, we have set the source terms Sa, S^^, and Sg to 
zero for simplicity but without loss of generality with re- 
spect to the main conclusions. If the shift is not outward 
pointing at the boundaries, and the initial data is chosen 
such that Ai = 0, Bi^ = 0, Ci = and Vkij = 0, then 



these results show that a solution to the first-order BSSN 
evolution equations will automatically satisfy these con- 
straints for all times. It follows that such a solution will 
also satisfy the original second order BSSN system. If 
the shift is outward pointing at a boundary, additional 
boundary conditions need to be specified in order to en- 
sure that these constraints propagate, see the discussion 
below Eq. (16). 
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On the other hand, numerical errors can trigger small 
violations of the constraints and these violations might 
grow in time. We can use the parameters k to help in- 
sure that the constraints are damped. As we show in the 
next subsection, k*^ = is required for strong hyperbol- 
icity, so let us consider k*^ = here as well. Now observe 
that with the standard gauge conditions the functions / 
and G are independent of 0. In this case the first two 
equations, Eqs. (18a)-(18b), decouple from the last two. 
With k" and sufficiently large, Ai and Bi^ should 
be damped. Next, observe that the fourth equation, 
Eq. (I8d), is independent of Ci. Assuming Ai and Bi^ 
are damped, the constraint 'Dijk should remain damped 
for sufficiently large constant . 

Finally, consider Eq. (f8c) with k"^ = 0. With the 
constraints Ai and Bi^ vanishing, this equation reduces 
to dtCi — CjjCi^ where Cp is the Lie derivative along the 
shift vector. It follows that the time evolution for Ci 
is simply a spatial diffeomorphism defined by the shift 
/3'. If initially the tensor components Ci are given small 
nonzero values due to numerical error, these errors should 
stay small as long as the spatial coordinates remain well 
behaved. 



C. Hyperbolicity 

The evolution equations (14)-(l-5) form a quasilinear 
first order system. 



dtu^ A{uy^^u + F[u), 



(19) 



F de- 



where the matrices A^, A^ and 

pend smoothly on the state vector u — 
{a,ai,l33,B^,l3i^,K, (/>, 0^ , 7^^- , Aij , %ij , A' ) . Such 
systems possess a local in time well-posed Cauchy 
problem if they are strongly hyperbolic, meaning that 
for each constant state vector u in an appropriate open 
neighborhood and each normalized covector there 
exists a symmetric, positive definite matrix H{u,n), de- 
pending smoothly on u and n^, such that H{u, n)A{uyni 
is symmetric. The motivation for this definition stems 
from the principle of frozen coefficients [65] in which 
the system (19) is first linearized about some smooth 
solution u, and then its coefficients are frozen at a 
specific point p of the spacetinie manifold. Denoting by 
u = u(p) the constant field that is obtained by freezing 
u at p, and by v the linearization of u, the system (19) is 
then replaced by the linear, constant coefficient problem 



dtv = A{uyd,v + T, 



(20) 



with T some constant vector. When T — 0, this sys- 
tem describes the evolution of small amplitude, high- 
frequency perturbations of the quasilinear system (19). 
Therefore, it is clear that a necessary condition for the 
well-posedness of the Cauchy problem for Eq. (19) is that 



the principal parts of all frozen coefficient problems, i.e. 
Eq. (20) with = 0, lead to well-posed Cauchy prob- 
lems. This turns out to be the case if and only if there 
exists a symmetric, positive definite matrix H{u,n) such 
that H{u,n)A{uyni is symmetric [65]. Provided H{u,n) 
depends smoothly on u and n, the principle of frozen co- 
efficients asserts that this is also a sufficient condition for 
the local in time well-posedness of the quasilinear prob- 
lem [65, 66]. 

The existence of the "symmetrizer" matrix H{u,n) im- 
plies, in particular, that the principal symbol A(u, n) :— 
A(uy"ni is diagonalizable and has a real spectrum for each 
u and n. Once this necessary condition has been verified, 
the symmetrizer H(u^ n) can be constructed by diago- 
nalizing A{u,n) = S{u,n)A{u,n)S{u,n)^^ with A(?2,n) 
a real, diagonal matrix, and then setting H{u,n) := 
{S{u,n)^^)'^ S{u,n)^^ . If S{u,n)^^ depends smoothly 
on u and n, this yields the required symmetrizer. The 
rows of S{il, n)~^u are the characteristic fields of the sys- 
tem (19), and the diagonal entries of A(m, n) are the cor- 
responding characteristic speeds. 

In our system (14)-(15) the principle part naturally 
splits into two blocks, one of them, the "gauge block" , 
comes from the evolution equations (14) for the 20 in- 
dependent variables a, a^, (5^, , /J^^, and K, and the 
other block, the "non-gauge block" , comes from the evo- 
lution equations (15) for the remaining variables. We 
first analyze the gauge block which is decoupled from 
the remaining block. Let us choose cr = 1. Through 
the replacements 9t 1— >■ ^ and di 1— >■ n.^, we find that the 
eigenvalue problem = A(u,n)v for this block reads 



ifJ- - /3n)a 
(/i - /3„)ai 
{fi - (5^)13 j 

(Ai - MB, 
(M - MK 



, 

—cP' JriiK + K^Uia , 
, 

4 
3 

a^GuiBj + K^UiP. 



4" 

3 ^ . 

'0 ' 



(21a) 
(21b) 
(21c) 

(21d) 

(21e) 
(21f) 



Here and in the following, the quantities d, /?', (j), jij re- 
fer to the frozen lapse, shift, conformal factor and phys- 
ical metric, respectively. Also, / = /(d, cj)) with similar 
definitions for G and H. We assume that /, G, and 
H are all positive. The covector rii is normalized such 
that j^^Uirij — 1. An index n refers to contraction with 
n* = 7'-'rij; for example, q;„ = n^ai. We have also used 
the frozen physical metric to lower indices: /3i = ^ijjS^ 
and = py^^kj- 

The characteristic fields and speeds for the gauge block 
are 
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AB 



ctA , Pa , Pau , a , f3ri , 



G 



(±) 



Pa 



■Vf 

4* (Aif 



3„ idVGTJ 
3„ ± aX . 



(22a) 
(22b) 

(22c) 

(22d) 



Here, indices A and B refer to contraction with unit vec- 



tors orthogonal to n% and we have set A = 
The characteristic fields are well-defined and indepen- 
dent from each other as long as 7^ /. This restriction 
on hyperbolicity, which is more explicity written as 



(23) 



is also required for strong hyperbolicity in the second 
order BSSN system [47]. 



The eigenvalue problem fiv 
gauge block is given by 



A{u, n)v for the non- 







= , 




- /3„)0» 


a 
~6 


(m- 


P 71)^ ij 


= , 




Pn)Aij 


a 
~2 


(m- 


Pn)lkij 


= -2c 


{^^- 


- °Pn)~^o 


— Pnj 



6 



rinij\ + e 



TF 



2ankAij + 2e '''^n^ [P(ij)] + K'^Ukjtj , 

'jn + T^rijPk' - UjK , 



(24a) 
(24b) 
(24c) 
(24d) 
(24e) 
(24f) 



r 



where A^. — e~^'^7fc^A^. This block consists of 32 inde- 
pendent variables: 1 variable </), 3 variables ^i, 5 variables 
7y- (since is the linearization of a metric with fixed 
determinant), 5 variables Aij (since Aij is symmetric and 
trace-free), 15 variables ^kij (since ^kij is symmetric and 



trace-free in i and j), and 3 variables Aj. 

The non-gauge block needs k"^ = to be diagonaliz- 
able. With = Q the characteristic fields and speeds 
are 



, Zq — (f>n — |A„ , (j)A , Jij 

y{±) ^ 7tf _ f^~tf 1 



Z,, = HA, - B,, 



lAij 



(AB) + 2 7„AB 



(±) _ 



y(.±) 
y nn 



nA 



= AnA - ^InA - ie~^^PAn T \ 



InnA 



2 'Inmi 



~ -f;OiA 
A„ — 20- 



M = 


/3n , 


(25a) 


M = 


/3„ ± d , 


(25b) 




/3„ ± d , 


(25c) 


M = 


/3„ ± d , 


(25d) 



where the superscript tf refers to the trace-free part in ^SabS'~''^ Acd- 



the transverse directions; for instance, A 



tf 

AB 



Aab ^ Provided the functions /, G and H depend smoothly 

on (q;,(/)) and satisfy the restriction (23), a smooth sym- 
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metrizer H(u, n) can be constructed from the characteris- 
tic fields as described at the beginning of this subsection. 

In the analysis above, we have assumed that all the 
algebraic constraints are identically satisfied, which is 
consistent with the evolution equations after the mod- 
ifications to Eq. (f5c) described in Section IIIB.^ If none 
of the algebraic constraints are enforced, then Eq. (24) 
yields {^^ - /3„)(f ^7*^) = 0, (m - PnW' A^) = and 

{fJ--^n){Y^lkij) = "fc[-2a(7*^Ay)-f KT(7*J7ij)] which is 
a weakly hyperbolic system. In order to obtain a strongly 
hyperbolic system one could enforce only the trace con- 
straint Ai^ = and replace Vkij by its trace-free part 
over ij in the right-hand side of Eq. (15c). In this case, 
7*''7ij and Y^lkij are characteristic fields with speeds 
/i = Pn and one has to perform the replacements 

Inn l-> g I 7n« - IpABO I 

and 

7„™^3(^W -^InAB^ j 

in the expression for Kin' in Eq. (25d). 



IV. NUMERICAL EXPERIMENTS 

Here we summarize results of numerical experiments 
of the first order BSSN formulation described in the 
previous sections, with different numerical approaches 
and codes, from more traditional ones for which there 
is more experience (finite differences with adaptive mesh 
refinement), to a promising approach that only recently 
is making its way into numerical relativity (Discontinu- 
ous Galerkin Finite Elements, restricted here to spherical 
symmetry) . 

In more detail, our approach and summary of numeri- 
cal experiments with the FOBSSN formulation is in the 
following order: 

1. Section IV A: Two Apples-with- Apples tests [67], 
as well as result from single and binary black hole 
moving puncture simulations using finite differ- 
ences with AMR. Most of our simulations show no 
signs of time or numerical instabilities. By time- 
stability in time dependent problems it is referred 
to the numerical solution not growing at any fixed 
resolution in time unless the exact solution does 
so. Numerical stability refers to the property that 
at any fixed time the errors in the numerical so- 
lution decrease with increasing resolution. In our 



^ These modifications do not change the principal part of the equa- 
tions when the algebraic constraints hold. 



simulations we have found the solution to be both 
time and numerically stable. 

The non-linear gauge wave test with a large am- 
plitudes shows a global time-instability (yet not 
a numerical one) that is expected; see [67]. This 
suggests that the addition of the extra constraints 
present when enlarging the system to a purely first 
order formulation does not trigger any obvious in- 
stability. The extracted gravitational waves are 
found to be consistent with simulations done us- 
ing the standard second-order BSSN formulation. 
In addition, the FOBSSN results are often more 
accurate than BSSN results using the same resolu- 
tion. 

2. Section IV B: With the standard BSSN gauge con- 
ditions a non-rotating black hole is driven to the 
trumpet solution [68]. However, this would require 
either the moving punctures technique [2, 3] or the 
turduckening one [56-58] . In the first case the equa- 
tions become singular at the puncture locations, 
which would be difficult to deal with using a very 
high order method such as those motivating the 
current paper. The turducken approach, on the 
other hand, smooths the solution inside the black 
hole while guaranteeing that the associated con- 
straint violations do not "leak" to the exterior of 
the black hole. As a first step in that direction we 
test a discontinuous Galerkin (dG) approach using 
the FOBSSN system for black holes in spherical 
symmetry, first using excision. 

3. Section IV C: As a final step, we perform turducken 
black hole dG simulations in spherical symmetry, 
both using FOBSSN and the standard second order 
formulation. We are able to perform long term and 
stable evolutions with the standard second BSSN 
formulation, but find numerical instabilities with 
FOBSSN. 

Discussions about all these experiments, their inter- 
pretation, and proposed next steps are discussed in Sec- 
tion V. Next we provide a somewhat detailed summary 
of these numerical experiments. 

A. Finite differences 

We have implemented the first-order system (14)-(15) 
using the Cactus framework [69, 70], and employing the 
Carpet adaptive mesh refinement (AMR) driver [71, 72]. 
We used the Mathematica package Krcoic [73, 74] to ex- 
pand the FOBSSN equations to C code, in the same man- 
ner as already for the McLachlan code [58, 75]. Both the 
Mathematica notebook as well as the resulting C code 
will be made available for public download as part of the 
Einstein Toolkit [76, 77] under the name Carlile. 

Our implementation supports arbitrary finite differenc- 
ing orders and time integration orders; below, we use 
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Code Name 


Symbol Eq. 


Value Comment 


harmonicN 


«V 


(14a) 


1 


(1 + log) 


harmonicF 




(14a) 


2.0 


(1 + log) 


ShiftGammaCoeff 


a^G 


(14c) 


0.75 


(std. choice) 


BetaDriver 


V 


(14d) 


1.0 




DAlphaDriver 




(14a) 


1.0 




DBetaDriver 




(14c) 


1.0 




DphiDriver 




(15a) 


0.0 


(not enforced) 


DgtDriver 




(15c) 


1.0 




sigma 


o 


(15f) 


1.0 





Robust Stability Test 



TABLE I: Constraint damping and related parameter settings 



fourth order accurate stencils and a fourth order Runge- 
Kutta time integrator. We use fifth order Kreiss-Oliger 
dissipation as well as fifth order spatial interpolation 
at AMR boundaries. We use buffer zones and tapered 
grids [58] to avoid time interpolation at mesh refinement 
boundaries. This makes all simulations fully fourth or- 
der convergent. The algebraic constraints 7*^Aij = 
and Y''lkij = are enforced every time the state vector 
is modified. However, 7 = 1 = 7 is not enforced, but is 
nevertheless assumed to hold throughout the implemen- 
tation. Our constraint damping and related parameter 
settings are listed in table I. We impose simple outgoing 
radiation (Sommerfeld) boundary conditions on all fields. 

a. Robust Stability Test. One of the most important 
and most fundamental test for a formulation of the Ein- 
stein equations and its numerical implementation is a ro- 
bust stability test, which can demonstrate linear stability. 
The simulation domain is initialised with Minkowski data 
plus a small amount of noise, and then let to evolve freely 
[67]. Here, we use a cubic domain with 40"^ grid points 
and periodic boundary conditions, and a noise amplitude 
of A = 10-6 in all BSSN or FOBSSN variables. Fig- 
ure 1 compares the performance of BSSN and FOBSSN, 
and finds very similar behaviour. In particular, the L2 
norm of the Hamiltonian constraint decreases steadily 
over time, indicating robust stability. 

b. Nonlinear Gauge Wave. A very demanding test 
is evolving a nonlinear gauge wave. This is a fully non- 
linear solution of the Einstein equations where the exact 
solution is known, as it is a flat spacetime in a complex, 
time-dependent coordinate system [67]. Here, we use a 
one-dimensional domain with 40 x 1 x 1 grid points with 
periodic boundaries, and evolve with the full 3D formu- 
lation. We test two large amplitude {A = 0.1) 
and a small amplitude {A = 0.01), employing the expsin 
form of the gauge wave. 

Figure 2 shows results from the large-amplitude case. 
This is a very demanding case that is known to go unsta- 
ble quickly for many formuations of the Einstein equa- 
tions [67]. Here we observe that both the evolutions 
with BSSN and the FOBSSN formulations break down; 
however, the FOBSSN evolution lasts for about twice as 
many crossing times. We also observe that the break- 
down mechanisms for BSSN and FOBSSN are different 



E 
ra 



0.01 



0.001 



0.0001 



1e-05 



E I I r 



~T 1 r 

BSSN 

FOBSSN 




100 200 300 400 500 600 700 800 
t[M] 

FIG. 1: Robust stability test comparing of BSSN and 
FOBSSN. A cubic domain is initialised with Minkowski data 
and a low level of noise in all variables, and then evolved with 
periodic boundary conditions. This tests linear stability of 
the formulation. The fact that the constraint violation de- 
creases indicates stability. Both BSSN and FOBSSN perform 
very similarly here. 



Non-Linear Gauge Wave, large amplitude (A=0.1) 



exact 
BSSN, t=70 
FOBSSN, t=1 20 




FIG. 2: Nonlinear gauge wave test, A = 0.1 (large amplitude), 
comparing of BSSN and FOBSSN. At t = 70 (35 crossing 
times), the BSSN solution has broken down (become irregu- 
lar) due to accumulation of numerical errors. The FOBSSN 
breaks down much later, shortly after t = 120 (60 crossing 
times); at t = 120, the FOBSSN solutions is still regular, and 
has only picked up a phase error and a global downwards drift 
in the metric. 



- the BSSN result develops high-frequency noise (de- 
picted), while in the FOBSSN result the metric drifts 
downwards, i.e. the proper size of the simulation domain 
decreases. 

Figure 3 shows results from the small-amplitude case. 
This is a less demanding case where most formulations 
of the Einstein equations can perform long-term evolu- 
tions [67]. After 100 crossing times, both the BSSN and 
FOBSSN result looks fine; however, the BSSN result ex- 
hibits a much larger upwards drift in the metric. 
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Non-Linear Gauge Wave, small amplitude (A=0.01 ) 



Total Mass of Puncture Black Hole 



exact 
BSSN, t=200 
FOBSSN, t=200 





BSSN 
FOBSSN 



30 40 50 60 70 
t[Ml 



80 



FIG. 3; Nonlinear gauge wave test, A — 0.01 (small ampli- 
tude), comparing of BSSN and FOBSSN. At t = 200 (100 
crossing times), both the BSSN and FOBSSN solutions are 
still fine. However, the BSSN solution has begin to drift up- 
wards much more than the FOBSSN solution. 



c. Single Puncture Black Hole. A much more in- 
teresting test of the FOBSSN formulation is evolving a 
puncture black hole. Here we choose a rotating puncture 
with total mass M — 1 and spin a = 0.7, set up via 
the TwoPunctures thorn [78]. These initial conditions 
are conformally flat and contain some gravitational radi- 
ation, and the black hole is expected to relax to a station- 
ary state after some time. In the figures below, we use 
a length unit M that corresponds approximately to the 
ADM mass of the system, which is Madm = 1.00252 Af. 
The black hole horizon has a coordinate radius of approx- 
imately 0.376 initially and 0.766 at late times. 

We employ eight levels of mesh refinement in a cu- 
bic domain, placing refinement boundaries at a; = 
[1, 2, 4, 8, 16, 64, 128] M, and placing the outer boundary 
at 258.048 M. The resolution on the finest level, which 
encompasses the horizon at all times, \s h — 0.032 M. 

Figure 4 shows the total mass of the black hole as 
calculated by the QuasiLocalMeasures thorn [79]. Af- 
ter an initial transient lasting about 20 M, the space- 
time becomes manifestly stationary. The angular mo- 
mentum (not shown) remains approximately constant at 
J = 0.701 ± 0.006 AP. Figure 5 shows a snapshot of the 
Hamiltonian constraint in this stationary state along the 
X axis si t — 76.8Af . As expected, the constraint viola- 
tion increases towards the black hole. Both BSSN and 
FOBSSN perform approximately the same except near 
the outer boundary, where FOBSSN seems superior. 

d. Inspiralling Binary Black Holes. As a more ad- 
vanced test, we also evolve inspiralling binary black holes, 
using the Rl configuration of [80, 81]. This configuration 
performs about 1.8 orbits prior to merger, with a common 
apparent horizon first found at roughly t = 160 M, where 
the ADM mass Madm = 0.966 M sets the scale. The 
initial individual black holes have masses Mi = M2 — 
0.505 M and have no spin. 



FIG. 4: Black hole total mass vs. time for a single punc- 
ture black hole with spin a = 0.7, comparing the accuracy 
of BSSN and FOBSSN. This is an initially non-stationary 
solution that evolves towards a trumpet solution. BSSN and 
FOBSSN perform very similarly here, and in particular, the 
moving puncture/turduckening approach to singularity han- 
dling seems to work fine for FOBSSN. 



Puncture Black Hole at t = 76.8M 
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FIG. 5: Hamiltonian constraint violation along the x axis at 
t = 76.8 M for a single puncture black hole, comparing the 
accuracy of BSSN and FOBSSN. The constraint violation 
increases towards the black hole (located at 2; = 0), where 
the horizon has a coordinate radius of about r — 0.766 M at 
this time. The constraint violation near x = 100 A/ is caused 
by outer boundary effects. FOBSSN seems to perform slightly 
better than BSSN in the bulk of the domain, and significantly 
better near the outer boundary. 



We use 9 levels of adaptive mesh refinement and placed 
the outer boundary at 320 M. The simulations were per- 
formed at two resolutions of h = Af/28.8 and Af = 
1/38.4, where h denotes the gridspacing on the finest 
grid, see also [80, 81] for comparison. We use fourth 
order accurate finite differencing stencils with lop-sided 
stencils for advection terms [82] and fourth order Runge 
Kutta time integration with Berger-Oliger sub-cycling in 
time. We do not employ tapered grids, using second or- 
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50 100 150 200 250 300 



time (M) 

FIG. 6: Comparison between simulations of the standard first 
order in time, second order in space implementation of BSSN 
and a fully first order reduction for a binary black hole inspiral 
system. The top panel shows the amplitude of the 1 = 2, 
m = 2 mode of ^'4 extracted at r = 50 at two resolutions 
h for both implementations. The bottom panel shows the 
difference between different resolutions. 

der time interpolation where necessary on mesh refine- 
ment boundaries. 

Due to the larger number of constraints in the first 
order formulation, one would expect a better accuracy 
in the second order formulation for the same number of 
grid points [83-86]. That seems indeed to be the case 
here: we find that the BSSN formulation allows us to 
use a lower resolution than the FOBSSN formulation to 
achieve time-stability. 

Figure 6 shows the amplitude of the £ = 2, m = 2 mode 
of the Weyl scalar ^4, extracted on a coordinate sphere 
with radius r = 50 M. The low-resolution FOBSSN sim- 
ulation is visibly different from the other simulations at 
the peak of the amplitude. However, the high resolution 
FOBSSN simulation agrees with both BSSN resolutions. 

Throughout the simulations, there is generally a good 
agreement between all runs. The lower panel of figure G 
shows the difference in the amplitude between high- and 
low-resolution runs for the BSSN and FOBSSN formula- 
tions, indicating that BSSN may have a smaller relative 
error. 



B. Discontinuous Galerkin 

Next we consider a dG scheme for the spherically re- 
duced first order BSSN system (see [51] for a dG im- 
plementation of the second order form of the BSSN 
equations). There are some important differences with 
FOBSSN in three dimensions, which arise when special- 
izing to spherical symmetry. First, the constraint A\ = Q 
is exactly satisfied by virtue of the spherically symmetric 
restriction. Second, terms proportional to a in Eq. (15f) 
are identically zero and so we set cr = 0. Finally, spherical 
symmetry is no longer associated with the obvious choice 

r^fe = and 7 = 1. As a consequence, using the fiducial 



covariant derivative will give rise to terms which feature 
7 and its derivatives. Our approach is to notice the fact 
that the covariant divergence I?i/3' only depends on 7 
and use the constraint 7 — 7 = to replace Dif5^ 
throughout BSSN system (4). Furthermore, we use r*" in 
place of A*", which results in a strongly hyperbolic spher- 
ically reduced BSSN system, explicitly given by Eq. (10) 
of [51]. For a complete discussion see Ref. [87]. 

We have discretized the first-order spherically reduced 
system with a nodal dG method [88, 89]. Similar to a 
multi-domain pseudo-spectral collocation method, a dG 
approach provides for a multi-domain treatment of the 
geometry where the numerical solution on each subdo- 
main is given by a (time-dependent) polynomial of arbi- 
trarily high order degree N . On every subdomain each 
component of the PDE is required to be satisfied in a 
suitable weak (integral) sense, yielding [N -\- 1) ordinary 
differential equations often known as Galerkin conditions. 
Adjacent subdomains are coupled in a stable manner 
through a suitable numerical flux term [88]. The result- 
ing scheme is nearly identical to the one presented in [51] 
with the notable exception of the absence of second or- 
der operators. Hence, we use the standard local Lax- 
Friedrichs form for the numerical flux [88], while in the 
second-order system, to which we will sometimes com- 
pare, a penalized central flux is used for a stable treat- 
ment second order operators (see page 13 of Ref. [51] for 
more details) . The integration in time is implemented us- 
ing the method-of-lines with a fourth order Runge-Kutta 
scheme. After each timestep an exponential fllter is ap- 
plied to the top two-thirds of the modal coefflcients to 
control alias driven instabilities. Furthermore, in our 
dG implementation of both the second and first order 
BSSN system we have empirically observed that the con- 
formal metric coefficients must not be filtered otherwise 
the scheme becomes unstable. A perhaps related obser- 
vation is that enforcing the constraint 7 = 7 triggers an 
instability at very early times. Neither this constraint 
nor the spherically symmetric version of Eq. (13) are en- 
forced in our dG implementation. 

All simulations presented next are for the Schwarz- 
schild metric in conformal ingoing Kerr-Schild coordi- 
nates'^ [51]. The source terms Sa, Sp, and Sg in the 
gauge equations (9) and (10) are chosen so that the nu- 
merical solution is time-independent. Typical choices for 
/, G, rj and a are used; in detail: / — 2/a, G = 3/4a~^, 
•q = 50, and terms proportional to a are identically zero. 
Furthermore, we set* H = e^'^jL and choose L 10 
such that the excision surface is not too close to r = 0, 
where field gradients are large. All damping parameters 



^ In these coordinates, at least initially, 7ij = diag(l, r^,r^ sin^ d) 
and so the algebraic constraint is 7 = 7 = r'* sin'^ 6. 

^ The dG code evolves the conformal factor x = e"'*'^. Neverthe- 
less, we continue to refer to the conformal factor as </> in this 
section. 
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H\\ at T = 3,000M 
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time 



FIG. 7; Discontinuous Galerkin evolutions of a black hole in 
spherical symmetry, using excision. See the text in Section 
IV B for more details. The last two resolutions have essen- 
tially reached double precision roundoff errors. 



k", , k'*' and are set to 20. We find that different 
values of L have negligible effect on the scheme's stabil- 
ity, and in particular no dependence on the location of 



the e 



4<j> 



2aL surface of potential weak hyperbolicity 



[cf. Eq. (23)]. 

The radial domain [0.4, 50]M is covered by 100 equally 
sized subdomains^ . We treat the inner boundary by exci- 
sion. At the outer physical boundary we specify the ana- 
lytic values for the incoming characteristic modes, which, 
for the spherically reduced system considered here, are 
given by Eqs. (17a-i) listed in [51] (Dirichlet conditions). 
Figs. 7 and 8 show that the scheme converges exponen- 
tially with N and is able to achieve very long run times. 
In an attempt to remove the slow growth in time for 
any fixed resolution seen in the Hamiltonian constraint 
we varied our numerical setup including the exponen- 
tial filter parameters (the number of filtered modes, the 
dissipation exponent, and which variables to filter), the 
timestep, the damping parameter rj, the numerical 
flux dissipation parameter, and the auxiliary field damp- 
ing parameters, without significant improvements. 



C. Turduckening 

Successful numerical evolution of binary black hole sys- 
tems require a suitable treatment of singularities. There 
are three distinct techniques used: moving-punctures [2, 
3], excision, and smoothing via turduckening [56-58]. 




This is far from optimal, since a better choice would be to have 
the size of the domains increase with radius. However, it suffices 
to make our point about stability and convergence. 



FIG. 8: Exponential convergence of Discontinuous Galerkin 
evolutions with polynomial order, see Section IV B. The norm 
used is the L2 one. 



State-of-the-art second order BSSN codes avoid the com- 
plications of excision, which require horizon tracking. A 
moving-puncture technique was used in our finite differ- 
ence implementation in Section IV A, while the dG code 
in Section IV C) each relied on excision. It has been 
shown that the usual gauge conditions are attractors of 
the trumpet solution [90-93] for which there is an in- 
coming characteristic mode even at the puncture [94] - 
generically we therefore do not expect an excision surface 
where no boundary conditions are required to exist. As 
the majority of BSSN implementations without excision 
have been thus far limited to finite difference methods, 
one wonders how other methods might deal with singu- 
larities. In this subsection we give a preliminary look at 
turduckening for the nodal dG code. 

We follow the turduckening technique described in 
Ref. [58] . Singular initial data in the interior of the black 
hole is replaced with smooth constraint violating data. 
The prescription for such smoothing used here is as fol- 
lows. If the computational domain is r e [0, R^g^^/M], 
we select a coordinate location rt inside the horizon and 
make the replacement r — f in the equations for the ini- 
tial data, where r is rigged to satisfy r(0) = tq, f{rt) = Tf , 
and f = r for r > rt. In addition, we require f to have a 
specified number of continuous derivatives (typically 8) , 
such that the turduckened and original data match to 
the specified degree of smoothness where they are joined. 
A polynomial f with these properties is constructed by 
solving a system of linear equations for the polynomial 
coefficients. 

In effect, our prescription stretches the physically cor- 
rect (non-singular) data for the region [tq , rt] over the 
turduckened region [0, rt]. This choice of initial data will 
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tern we have two distinct choices for calculating the aux- 
ihary variables in the turduckened region of the initial 
data: calculating the analytic derivatives of the fields 
at the turduckened grid points or applying the numer- 
ical derivative operator to the turduckened fields. In 
the first case the auxiliary constraints are violated since 
the auxiliary fields correspond to derivatives of the non- 
turduckened fields. In the second case the auxiliary con- 
straints are satisfied but the turduckened initial data no 
longer represents the physically correct data for [ro,rt] 
stretched over region [0,rt]. We experimented with both 
choices and found that the region of constraint violation 
is not guaranteed to shrink when using the first choice 
in which the auxiliary constraints are violated. Fig. 10 
documents a typical comparison with turduckening pa- 
rameters rt = .3M and vq = AM. We note that our 
observations are influenced, for example, by the source 
terms Sa, S^, and Sg. 



FIG. 9: Evolutions of a black hole in spherical symmetry using 
the standard second order in space BSSN formulation and a 
dG scheme with turduckening, see Section IV C for details. 
The norm used is the L2 one. 

naturally be constraint violating. As the constraint sys- 
tem's wavespeeds (see Section IllC and Ref. [58]) are 
not superluminal, these violations remain "trapped" in- 
side the horizon for all future times. Furthermore, for the 
second order BSSN system, Ref. [58] found that the re- 
gion of constraint violation quickly shrinks relative to the 
numerical grid. We experimented with turduckening the 
second order dG scheme described in [51], and found that 
the region of constraint violation quickly shrinks with 
this scheme as well. These simulations used a grid with a 
larger outer boundary and staggered domain sizes: 1 sub- 
domain [0,rt = A]M comprising the turduckened region 
with To = .IM, 3 subdomains in [.4, 1.5]M, 6 subdomains 
in [1.5,10]M, and 12 subdomains in [10, 100]M. Other- 
wise the same numerical settings as in IV B, although 
here we use frozen outer boundary conditions on all fields 
and gauge source terms chosen so that Eqs. (9, 10) are 
initially time-independent. Fig. 9 shows the scheme is 
stable and whenever t > M converges towards a time- 
independent solution exponentially with N. 

Furthermore, the technique is observed to be robust 
for a variety of numerical parameter choices and domain 
decompositions. Results from our second order BSSN 
dG code suggest the turduckening technique to be appli- 
cable beyond finite difference schemes. Nevertheless, we 
have been largely unsuccessful at achieving robust sta- 
bility turducken tests for FOBSSN. A typical evolution 
lasts on the order of tens to hundreds of M, although 
some low resolution runs can last into the thousands of 
M before crashing. 

The presence of extra auxiliary constraints (12d)-(12a) 
presents a genuine difference between turduckening a first 
and second order BSSN system. In the first order sys- 



V. COMMENTS 

The goal of this paper has been a first step towards 
combining the robustness and simplicity of evolutions of 
the BSSN formulation of Einstein's equations, most no- 
tably being able to avoid the complications of excision, 
with very high accuracy numerical schemes - those being 
a multi-domain pseudo spectral collocation method and 
a discontinuous Galerkin method. Furthermore, any of 
these approaches would allow, due to their memory ef- 
ficiency and the speedup of GPUs (Graphics Processing 
Units), to run binary black hole simulations on a single 
GPU, thereby avoiding the bottleneck of PCIe commu- 
nication between CPUs; see for example [95]. 

For this purpose we derived and analyzed the hyper- 
bolicity, characteristic variables and constraint propaga- 
tion of a fully first order BSSN formulation of the Ein- 
stein equations with optional constraint damping terms, 
FOBSSN. Unfortunately, we have not been able to derive 
a symmetric hyperbolic formulation, but only a strongly 
hyperbolic one. It is known that in more than one spa- 
tial dimension, strong hyperbolicity, even with maxi- 
mal dissipative boundary conditions, does not guaran- 
tee well posedness of the initial-boundary value prob- 
lem [96]. Yet, in our numerical experiments we have been 
able to carry out binary black hole simulations using our 
FOBSSN system, finite differences and adaptive mesh re- 
finement, without any need for fine-tuning and no obvi- 
ous signs of time instability (convergent errors that grow 
in time) or numerical instability (errors that get larger 
at higher resolutions at any fixed time). Most notably, 
the presence of the extra constraints in FOBSSN seem to 
cause no problems. 

Next natural steps would be three-dimensional discon- 
tinuous Galerkin evolutions of FOBSSN using, for ex- 
ample. Hedge [97], and implementation of more sophisti- 
cated boundary conditions, such as those of Ref. [Gl]. 
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\-H\\ at T =0M: 



10' 



10 



10 



10 



Auxiliary constraints satisfied 
Auxiliary constraints violated 




|7-i:|| at T =im: 




Auxiliary constraints satisfied 
Auxiliary constraints violated 



FIG. 10: Snapshots of the Hamiltonian constraint for turduckened initial data which satisfies (black line) and violates (dashed 
green line) the auxiliary constraint conditions (12d)-(12a) which arise in the first order reduction. Notice that only in the 
first case the constraint violating region "shrinks" in time, see Section IV C for more details. The solid vertical line marks the 
location of the event horizon, and the dashed one the location of the outermost radius with purely outgoing modes (where 
excision could be performed). The large red asterisks on the horizontal axis mark the location of each subdomain boundary. 
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Appendix A: The covariant BSSN system 

In this appendix we sketch the derivation of the covari- 
ant BSSN system, Eqs. (4). The derivation follows from 
the analysis in Ref. [59] by setting the determinant of the 
conformal metric to 7, and choosing the trace of Aij to 
vanish. For simplicity we assume that the fiducial fields 7 

and r jk are constructed from a time-independent met- 
ric 7jj. Unlike in the main body of the paper, here we 
do not assume that the fiducial metric is flat. 

Begin with the evolution equations for the physical 
spatial metric and extrinsic curvature, 

5j_7»j = ~2aKij , (Ala) 

-D.DjU , (Alb) 

where i?^ and Di are the Ricci tensor and covariant 
derivative for 7^. Now let d± = dt — Cp act on the 
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BSSN variables 4> and 7^, which are defined in Eqs. (1). 
The right-hand sides of these equations are written in 
terms of BSSN variables by inverting the definitions (1): 

7„- = e4*7,, , (A2a) 
K,, = e^^i,, + ^j,jK . (A2b) 

The results are identical to Eqs. (4a) and (4b), respec- 
tively. 

The derivation of the equation of motion (4d) for Aij 
follows the same pattern; apply d± to Aij in Eq. (Id), use 
the evolution Eqs. (Al), then replace the physical metric 
and extrinsic curvature with the BSSN variables through 
Eqs. (A2). In this case we must also write the physical 
Ricci tensor Rij in terms of conformal variables. Insert 
the relation 

r'jfc - r^k + 2(<5pfc0 + SID,^ - %uD'<l>) (A3) 

for the ChristofFel symbols into the definition of the Ricci 
tensor. This yields the splitting (5) between the confor- 



With the 



If the fiducial metric is flat, as assumed in the main body 
of the paper, then the fiducial Riemann tensor term on 
the right-hand side vanishes. The result (6) is obtained 
by replacing AF* with the new variable A* and dropping 
the fiducial Riemann tensor. 

To obtain the equation of motion (4c) for we first 
let d± act on if = -f^^Kij, using the results (Al). The 
right-hand side is simplified by adding —aH, where T-L = 
— KijK^^ +R is the Hamiltonian constraint. Equation 
(4c) then follows after using the inverse relations (A2) to 
write the result in terms of BSSN variables. 

The conformal connection vector is defined in Eq. (3). 
To derive the evolution equation (4e) for A% we first let 
the operator d±_ act on AF' = j^'^AT^jk with AT^jk ex- 



mal Ricci tensor Rij and the terms (7) that depend on 
the conformal factor (j). 

The derivation of the identity ((>) used for the confor- 
mal Ricci tensor is somewhat tedious. Beginning with 
the definition 

Rij = dkT^ij — diV^ jk + T^ijT\i — T^iiT^ jk (A4) 

it is staightforward to show that the difference between 
the conformal and fiducial Ricci tensors is 

Rij — Rij — DkAT^ij — DiAr'^jk 

+Ar\,Af 'ki - Af\iAr',k ■ (A5) 

One can also show that 

Af = (Dnki + Dk^ji ~ Dajk) , (A6) 

and derive the useful relations Dk'jij — 2AF(ij)fc and 
75fc7^-'' = -2Af(^-?)fe. With these results, the first two 
terms in the difference (A5) become 



(A7) 



(A8) 
I 

pressed as in Eq. (A6). This generates several terms of 
the form d±(Di'^jk)- Using the fact that Lie derivatives 
and partial derivatives commute [98] one can write these 
as 

d±(Dajk) = D^id^jjk) + 2{C,3T\(^^)%i . (A9) 
Now use the identity [99] 

CpT'jk = DjDkP' - ItjkiP' , (AlO) 

which is straightforward to verify. The result of this cal- 
culation for 9_lAA* is an expression in which the opera- 
tor acts only on the conformal metric 7ij. Using the 
equation of motion (4b), we find 



DkAf^,, ~ AAP.fc = -\i'''DkDi% + D(,(Af,)fe'=) - {D(,^'''){Dkij)i) - \{Dki'''){Di%,) 

+ {Dkf )(D^af)i) - R,j - f'%,i Ji,)kr ■ 



definition Af * = 7^''Af "j^, the conformal Ricci tensor from Eq. (A5) becomes 

.^j,*;'7i_ 7i_;r, 1 ~ri AJ^k zMz, 15 rn , ~,kl 

I ' 

-7"(2Af™fc(,Af,)„, + At"\kAf^ji) 



R^3 = ^^f'DkDa^j + Jki^D^^AT" ~ f'jrn^,Rf)kr + 7"' AF^M AF(,,), 



9i(Ar) = ^'"D.Dkfi' - f'R^kiP' - -^D, (av^i'^) + r^j'" AV ,kDif3' + ^D^Okp") . (All) 



Next we add the term 2aA^^ which is proportional to the momentum constraint (8b), to obtain 
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2A'''dka + 2aA"ArV, + l2aA'''dk(t) - -aD'K 



(A12) 



I 

Now use the definition (3) to replace d±AT^ with d±A^. Appendix B: Fundamental variables in terms of 

The result, assuming the fiducial metric is flat so that characteristic variables 

R^jki vanishes, is the equation of motion (4e). 

The first-order BSSN variables in the gauge block can 
be obtained from the characteristic variables using the 
formulas 



Ba 

PnA 

K 

B„ 



(+) 



G 



(-) 



Pa , 



^A ^ 



(-) 



2 V / Ai^ 



PG 



AH 
3(A2 - /) a 



1 



-PabS^'' + 



\G 



2 V ) 3(A2 - f)\ a^X J 



(Bla) 

(Bib) 
(Blc) 

(Bid) 
(Ble) 

(Blf) 



We then apply the relations 



Ui = UiUn + 7i aA , 
Pij — PnnTliT^j + T^'i'^j^ PnA 

+%^njPA7i + li^lj^ Pab 



(B2a) 



with similar expressions for Pi and Bi. 



(B2b) For the non-gauge block, the inverse transformation is 
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A% = \ [V^/^ + Vjcj) + g + I e-'^PlL^ , (B3a) 
iJab = - (V\V Vjcj) , (B3b) 



Aa = ^{Za + Ba) , (B3c) 
H 

AnB = I (V!:P + Vy^^) + ^ 7ni3 + 4 e'^'' Psn , (B3d) 



InnB = - - + e"""^ (^B " 20i3 - 4^5 j , (B3e) 

A„ = 1(Z„ + S„) , (B3f) 
H 

0n = Zo + ^ A„ , (B3g) 

Ann = ^ + KL")) + g 7«n + ^ - \ e'^^/J^B^^^ , (B3h) 

7n„„ = - - KL;') + \ (a„ - 20„) . (B3i) 
The tensor components Aij and 7^^^ are then reconstructed as 

Aij = Ann Q - ^7y^ + ^i^o^AnB + + ^^^j^Aab , (B4a) 

7.. = %nnn, [l n.^n, - i^.) + n,n.^,-^nnB + n,n,7.-7™B + ^,7/7,^:^.. + 7.^7^. • (B4b) 
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